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ABSTRACT 

The majority of stars form in star clusters and many are thought to have planetary 
companions. We demonstrate that multi-planet systems are prone to instabilities as a 
result of frequent stellar encounters in these star clusters much more than single-planet 
systems. The cumulative effect of close and distant encounters on these planetary 
systems are investigated using Monte Carlo scattering experiments. We consider two 
types of planetary configurations orbiting Sun-like stars: (i) five Jupiter-mass planets 
in the semi-major axis range 1 — 42 AU orbiting a Solar mass star, with orbits that are 
initially co-planar, circular, and separated by 10 mutual Hill radii, and (ii) the four gas 
giants of our Solar system. We find that in the equal-mass planet model, 70% of the 
planets with initial semi-major axes a > 40 AU are either ejected or have collided with 
the central star or another planet within the lifetime of a typical cluster, and that more 
than 50% of all planets with a < 10 AU remain bound to the system. Planets with 
short orbital periods are not directly affected by encountering stars. However, secular 
evolution of perturbed systems may result in the ejection of the innermost planets or in 
physical collisions of the innermost planets with the host star, up to many thousands 
of years after a stellar encounter. The simulations of the Solar system-like systems 
indicate that Saturn, Uranus and Neptune are affected by both direct interactions 
with encountering stars, as well as planet-planet scattering. Jupiter, on the other 
hand, is almost only affected by direct encounters with neighbouring stars, as its mass 
is too large to be substantially perturbed by the other three planets. Our results 
indicate that stellar encounters can account for the apparent scarcity of exoplanets 
in star clusters, not only for those on wide-orbit that are directly affected by stellar 
encounters, but also planets close to the star which can disappear long after a stellar 
encounter has perturbed the planetary system. 

Key words: planets and satellites: dynamical evolution and stability — planetary 
systems — open cluster and associations: general 



1 INTRODUCTION 

Since the disco v ery of the first exoplan ets by 

IWolszczan fc Frail (jl992h and iMavor fc Queloj (|l995h . 
hundreds of new exoplanets have been discovered, and 
thousands of exoplanet candidates have been identifie d 
using the Kepler observatory (e.g., iBatalha et all l201ot ). 
A substantial fraction of the kn own planets ar e part 
of a multi-planet system (e.g., ICumming et all 120081 : 
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iBorucki et ail l201ll : iLissauer et ail 120121 ). iFabrvckv et al.l 
i|2012l ) indicate that a large fraction of these exoplanet 
candidates are genuine, and that the majority of the planets 
in a multi-planet sy stem orbit their host star in nearly the 
same plane (see also iFigueira et aill2012l ). 

Most stars, and therefore most pl anetary systems, are 
thought to form i n star clusters (e.g.. lLada fc Ladall2003l ; 
IClarke et al1l2000l ). a lthough some may f orm in relatively 
isolated regions (e.g., iBressert et all 120101 1 . A certain frac- 
tion of these star clusters may survi ve for a long time, but 
many dissolve within 20-50 Myr (see lde Griis fc Parmentierl 
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120071 : IPortegies Zwart et al.l [2010L and references therein). 
Several studies have suggested that even o wn Solar system 
may h ave formed in an op e n clu ster (e.g.. IPortegies Zward 
2009). Dukes fc Krumholj ||2012h . on the other hand, ar- 
gue that current observations do not rule out that the Solar 
system may have fo rmed in a mas sive, but short-lived, star 
cluster. We refer to lAdams! (|2010h for an extensive review 
on the birth environment of the Solar System. In such dense 
stellar environments, dynamical encounters between stars 
are frequent and can perturb the orbit of binary stars and 
planetary system significantly. 

As star clusters are the likely birth places of the 
majority of the planet-hosting stars in the Solar neighbour- 
hood, numerous exoplanet surveys have been carried out in 
open clusters and globul ar clusters, but th eir results have 
been disappointing (e.g.. [Zhou et all |2012| . and references 
therein). Very few or no close-in planets with orbital period 
range of 1 — 16 days have been dete c ted in 47 Tucanae 
and u> Centauri dWeldrake et all 120051. 120081), the Hyade s 



i|Guenther et al.ll2005l). NGC 6397 jNascimbeni et~ai]|2012l) 
NGC 6791 ilMontalto et al l 120071) . M37 (iHartman et all 
I2009T). NGC 679 1 dMocheiska et all |2005|), NGC 2158 
ilMocheiska et all I2006D. NGC 7 789 |Bramich fc Home] 
l2006h. NGC 1245 dBurke et all I2006T ), and NGC 7068 
I Rosvick fc Robbl I 20061 ). However, a number of p lanet s 
hav e been found in open clusters: lLovis fc Mavorl (|2007r ) 
and ISato et alj (120071 ) both repo rt a planet in a wid e orbit 
around a giant star. In addition, iQuinn et al.l l|2012l ) report 
the discovery of two hot Jupiters in Praesepe, making 
these the first hot Jupiters to be detected orbiting a main 
sequence star in a star cluster. Finally, pr otoplanetary disks 
have been detected in star clusters (e.g., lHillenbrandll2005l . 
and references therein), hinting at planet formation. These 
observations demonstrate that planets can form in star 
clusters, but it leaves the question why so few have been 
discovered. 

The paucity of planets in star clusters may partially 
be explained by the suppression of planet formation in 
star clusters. Jupiter-mass planets are less likely to form in 
the neighbourhood of massive stars in star clusters due to 
the evaporation o f circumstellar disks by photo-evap o ration 



ilArmitagel 120001; iHoldert et all l201lh. lOlczak et all 
l2008h . iForgan fc Ricd (|2009l ) and iLestrade et all 



showed that stellar encounter s reduce planet formation due 
to disk destruction, although iThies et all (|2005l . l20ld ) find 
that encounters between circumstellar disks and other stars 
may actually enhance planet formation. 

After the formation process, a newly-formed planetary 
system and its cometary population remain prone to dis- 



cluster (e.g., Adair 


is & Laughlinl 200ll: lAdams et al.ll2006al: 


lAdams & Laughlin 


2006; Brasscr ct al. 201§). The impor- 



tance of the encounters is primarily determined by the lo- 
cal stellar density, the collisional cross section of the plan- 
etary system, and the timescale at which the system is ex- 
posed to externa l pertu rbations. iBonnell et all l|200ll ) and 
ISmith fc BonneDI l)200ll ) show that planetary systems similar 
to our own Solar Systems are li kely to survive i n ope n clus- 
ters, but not in globular clusters. iFregeau" et al.l f2006) study 
star-planet systems in star cluster environments and derive 
dynamical cross-sections and the location of the hard-soft 
boundary for star-planet systems, showing that these sys- 



tems are more prone to disruption than previously thought. 
Flybys can cause immediate ejections, but can also trigger 
instabilities in multi-planet systems that can lead to the 
ejection of one or more planets millions of yea rs after the 
encounter occurred (|Malmberg et alj|2007l . l201ll ) . A system- 
atic study on the evolution of s ingle-planet sys t ems i n glob- 
ular clusters was carried out bv lSpurzem et all l|2009l ). They 
find that close-in planets are difficult to disrupt, but their 
eccentricities are excited by weak encounters, which may 
ultimately result in their migration as a result of tidal evo- 
lution. 

Given that most star clusters may initially have a cool 
virial state and a fractal structure, close encounte rs may be 
more frequent during the ear ly phase of relaxation (jWoolfsonl 
120041 : IParker fc Quanzll2012l ). A highly inclined companion 
could excite plane ts to high ec centricity orbits through the 
Koza i mechanism ( Kozai 1962) as a result of secular interac- 
tion (|Wu fc Murravll2003l ). Studies on the secular instability 
of several isolated multi-planet systems caused by angular 
momentum deficit (AMD) and mean motion resonance show 
that these non-circular and non-co-planar orbits generated 
by the perturbation with random resonances would be un- 
stable (|Laskar fc Correia II2OO9I : lLaskar Il2000l ). 

The evolution of compact multi-planet systems in star 
clusters, however, remains poorly understood, although pre- 
vious studies have already indicated that the effect of planet- 
planet interaction s may be at least as important as the effect 
of st e llar flybys (lAdams et al.l l2006al; lAdams fc La ughlin 



120061 ; iMalmberg et al.ll201ll ; iBolev et al.1 12012| ). It is diffi- 
cult to model the evolution of these in star clusters using 
direct Af-body simulations, because of the enormous range 
of temporal and spatial scales between individual planetary 
systems and their host star clusters. In addition, the number 
of weak stellar encounters is large and numerical err ors in 
planetary systems are cumulative (e.g., Aarseth 2003). Pre- 
vious studies have therefore often focused on single-planet 
systems, which can be modelled as perturbed two-body sys- 
tems, or on planetary systems in low-mass, short-lived star 
clusters. 

In this article we present Monte Carlo scattering ex- 
periments to study the evolution of multi-planet systems in 
dense stellar environments, modelling all encounters within 
1000 AU during the time a typical planetary system spends 
in a dense open cluster environment. We find that the effect 
of flybys is important for the evolution of a multi-planet 
system. Interestingly, we demonstrate that the evolution of 
short-period planets is substantially affected by the pres- 
ence of outer planets. The article is organised as follows: we 
describe the methods and assumptions in § [5] The results 
and analysis are presented in § [3] and finally we draw the 
conclusions in § [4] and discuss our results in § [5] 



2 METHOD AND INITIAL CONDITIONS 

Under some circumstances it is possible to run full A-body 
simulations of dense star clusters with planetary systems. 
For single-planet systems, this is relatively easy, since two- 
body systems can be regularised. Spur zem et all (|2009l ). for 
example, used this method for a comprehensive study on 
the dynamical evolution of single-planet systems in popu- 
lous star clusters. The situation is substantially more com- 
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plex for multi-planet systems. Individual planetary motions 
require integration using time steps which are many orders 
of magnitude smaller than those used for the stellar dynam- 
ics. In addition, the integration errors in planetary systems 
are cumulative. As the secular behaviour of a multi-planet 
system is of utmost importance to its dynamical fate, it is 
currently impractical to use full iV-body simulations. Note 
that several attempts have been made to carry out such 
simulations with relati vely low-mass star cl usters and on 
shorter timescales fe.g. jMarmberg et al . 2011). To study the 
evolution of multi-planet systems in dense clusters on long 
timescales, however, we need to resort to alternative tech- 
niques. 

We use a Monte-Carlo scattering approach to simulate 
the environment in which most of the planetary systems are 
born. We model the evolution of the newly generated plan- 
etary systems at high accuracy and attempt to reproduce 
their environments more realistically than in previous stud- 
ies. Instead of carrying out time-consuming, and for our pur- 
pose inaccurate, direct TV-body simulations of star clusters 
with planetary systems, we focus only on individual plan- 
etary systems and their encountering stars, while ignoring 
the more or less homogeneous effect of other surrounding 
stars in the cluster. 

We first calculate the encounter rate for a typical star in 
an open cluster environment. Based on the total number of 
stars that approach the planetary system within a distance 
of 1000 AU during its lifetime in the star cluster, we let the 
encountering stars approach the planetary system one by 
one at a Poissonian time interval, derived from the average 
encounter frequency. 



2.1 Modelling the encountering stars 

We consider the evolution of planetary systems in open clus- 
ters similar to the Orion Nebula Cluster (ONC), which has 
roughly N = 2 800 stars wi th a one dimensional velocity dis - 
persion a = 2.34 kms" 1 ijflillenbrand fe Hartmann Nl998l ) 
and a half-mass radius of i?hm = 0.5 pc. Within the half- 
mass radius, the average distance between two stars is ap- 
proximately RhmN^ 1 ^ 3 . We assume that the star cluster 
is in virial equilibrium and ignore the effects of mass seg- 
regation and evaporation. The lifetime of an open cluster 
depends on the properties of the star cluster and its envi- 
ronment. Open clusters disso lve within 10 7 — 10 9 years (e.g., 
Ide Griis fc Parmentier]|2007l ) and we therefore adopt a total 
integration time of 100 Myr. 

The properties of the encountering stars are modelled 
by drawing (i) the relative velocity at infinity Voo, (ii) the 
impact parameter b, (iii) the mass M e of the encountering 
star, and (iv) the encounter rate V from the distributions 
that are representative for open clusters. For simplicity, we 
assume in this study that these four parameter distributions 
are mutually independent, and independent of time. All en- 
countering stars are initialised to approach from random 
directions. Our choices for these four parameters, as well 
as several other environmental quantities are summarised in 
Table E 

The mas s es Me of the encountering stars are drawn from 
the lChabrier] {2003) normal mass distribution, 



/(logM e )oc 



exp 



(1) 



where /i — log(0.2) and <r M = 0.55 in the mass range 
0.08 — 5Mq. Although more massive stars can cause sub- 
stantial damage to planetary systems during their encoun- 
ters, we ignore them since they are relatively rare and short- 
lived. As mentioned above, we assume that the mass of the 
encountering stars is uncorrelated with their impact param- 
eter, their velocity at infinity, and time. In other words, we 
ignore the effects of mass segregation and stellar evolution. 

The encounter velocities are derived from the veloc- 
ity dispersion of the star cluster. For the iPlummerl (Il91ll ) 
model, its relation to the total mass M c i us t C r and half-mass 
radius Rhm of the cluster is given by 



GM a 



r/R h 



(2) 



where r\ ~ 9 . 75, and G is the gra vitational constant (e.g., 
ISpitzeij [19871 ; IHeggie fc Hutll2003l ). In our simulations we 
are interested in the distribution of velocity differences Doo 
between two stars, which can be obtained from the Maxwell- 
Boltzmann distribution, 



f(Voo) = 



20Fa : 



■ exp 



~4cr 2 



(3) 



where the mean velocity difference between two stars is 
(voo) = 2<7\/27i- -1 = 3.74 kms -1 , and a is the one- 
dimensional velocity dispersion. As all planets have masses 
much smaller than either of the two stars involved, the 
stars follow near-hyperbolic trajectories. The relative veloc- 
ity v(r) between the two stars as a function of their separa- 
tion r(t) can therefore be approximated with 

i ^ 2GMt , ^ 
r(t) 

where Mt = M c + M e denotes the combined mass the bod- 
ies involved. The hyperbolic orbit can be described using 
its impact parameter b and its velocity at infinity Vaa. The 
impact parameter b is drawn from the distribution 

26 



Mb) 



(5) 



We only consider encounters with impact parameter < b ^ 
frmax, where 6 ma x is the maximum impact parameter, and we 
ignore the effect of stars with an impact parameter larger 
than brriaY. lAdams fc Laughlin] |200ll . [20061 ) have shown that 
the most important parameter that determines the effect of 
a stellar encounter on a planetary system is their distance 
of closest approach p. Its relation to the impact parameter 
b can be expressed as 



1 + 



1GM T 



(6) 



Typical encounters have (Mt) ~ 1.55M© and (voo) = 
3.74 kms - . For these interactions, b > 1094 AU corre- 
sponds to a distance of closest approach of p > 1000 AU. In 
our simulations we ignore the weak perturbations induced 
by stars with p > 1000 AU. 

Now that we have drawn Voo, b, and M e from their re- 
spective distributions, we can calculate the hyperbolic semi- 
major axis 
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Quantity 


Value 


Number of stars N 


2800 


Initial mass function 


Chabrier 0.08 - 5M 


Half-mass radius -R^m 


0.5 pc 


Velocity dispersion a 


2.34 kms" 1 


Closest approach p (AU) 


- 1000 AU 


Impact parameter b (AU) 


f b (b) oc b; fe max = 1094 AU 


Mean velocity difference »oo 


3.74 kms" 1 


Average encounter interval 


(t) = 1.7 Myr; Poissonian 


Total integration time 


T = 100 Myr 


Encounter number 


58 


Impact direction 


Random 



Table 1. Initial properties of the example open cluster and stellar 
encounters. 



GM T 



i4 Ve 2 - 1 e - 1 
and the hyperbolic eccentricity 



e=\ 1 + 



( bvl 



\GM T 



(7) 



(8) 



where p is the distance of closest approach. Finally, the ve- 
locity at periastron v p is given by 

2 GMWe + lX 

The encounter rate F is calculated using T = no A, 
where n is the local stellar density, a the velocity disper- 
sion and A the cross section of the encounter, which can be 
approximated with 

„.2 



A = ntf 



irp 



1 + 



= np I 1 + 



e+2 
1 



(10) 



The time r between two subsequent interactions is drawn 
from a Poisson distribution 



P(N(t) =r) = 



(tr- 



(ii) 



where (r) = is the average time between two subsequent 
encounters. For the star cluster listed in Table Q] we obtain 
(r) w 1.7 x 10 year, corresponding to a total number of 
58 stars that approach the planetary system within p — 
1000 AU. 



2.2 Initial conditions for planetary systems 

We model two sets of planetary systems, which we refer to 
as model 1 and model 2, respectively. The initial conditions 
for the planetary systems are summarised in Table [21 We 
adopt these two sets of planetary configurations in order 
to obtain a better understanding of the effect of different 
configurations of multi-planet systems. 

For model 1 we follow the EMS (equal- mass, equal- 
relativ e-separation systems) prescriptions of IZhou et al.l 
i|2007l) and we adopt a central star mass of IMq. Planet 
formation simulations have indicated that when the mass of 
a gas giant reaches around a Jupiter mass (lMj), its oli- 
garchic gr owth terminates and the planet obta ins its isola- 
tion mass l|Kokubo fc 7dall2002l ; llda fe Linll2004 ) . Therefore, 
we assign all planets a mass of lMj. These simulations also 



Model 


Quantity 


Value 


1,2 


Central star mass 


M c = IMq 


1 


Number of planets 


N p = 5 


1 


Planet mass 


M p = 1 Mj 


1 


Minimum semi-major axis 


Cmin = 1 AU 


1 


Mutual separation 


fc = 10 


1 


Eccentricity 


e = 


1 


Inclination 


i = 0° 


1 


Initial orbital phase 


0° < <f> < 360° 


2 


Number of planets 


N p = 4 


2 


Orbital parameters 


Jupiter, Saturn, 






Uranus, and Neptune 



Table 2. Initial conditions for the planetary systems used in our 
simulations, for model 1 and model 2, respectively. Note that 
model 2 represents our Sun with the four gas giants in our Solar 
system. 



indicate that most planetary embryos are initially separated 
by 10 — 12 mutual Hill radii (Rh) after they have depleted 
the planetesimal populations in their neighbourhoods. We 
therefore assign our planets orbits which are initially sepa- 
rated by k = lOiiff. The mutual Hill radius for two neigh- 
bouring orbits with semi-major axes ai and 02 is defined 
as 



Rh — 



01 +a 2 /Mi +M 2 



/ Mi + A/ 2 \ 
V 3M C J 



1/3 



(12) 



where Mi and M 2 are the masses of the two planets and M c 
is the mass of the central star. 

For the innermost planet we adopt a semi-major axis of 
flmin = 1 AU, which is within the ice line for solar-type stars, 
and accounts for a modest amount of migration. The semi- 
major axis of the subsequent planets is determined by the 
value of k. In our case we choose N p — 5 planets, resulting in 
semi-major axes for the planets of a — 1, 2.6, 6.5, 16.6, and 
42.3 AU, where the outermost semi-major axis reaches the 
limit at which planetary systems are able to form planets 
through the core-accretion process. All planets are assigned 
zero initial eccentricities (e = 0) and inclinations (i = 0°), 
and the initial orbital phase of each planet is drawn ran- 
domly from (j) e [0°,360°]. 

For model 2 we adopt a Solar-mass star with four gas 
giants: Jupiter, Saturn, Uranus, and Neptune. Each of these 
planets is assigned their mass and present-day orbital ele- 
ments identical to their counterparts in our Solar System. 
Note that the encountering stars come from all directions, 
so the inclinations of all orbits will gradually increase over 
time. 



2.3 Simulations 

The simulations are carried out by combining the strengths 
of two integration packages. During an encounter with the 
planetary system, the package CHAIN is used, while the 
MERCURY package is used to model the long-term secu- 
lar evolution of the system in be tween two subsequent en- 
counters The CHAIN package l|Mikkola fc Aarsethl Il993l ; 
lAarseth|[2003l) uses chain regularisation and is particularly 
useful for the treatment of energetic multi-body encounters. 



The integration is carried out quickly, and energy is con- 
served at a level of \AE/E\ < 10" 13 . Although CHAIN can 
very accurately model short-term energetic encounters in 
small- iV systems, it is not optimised for long-term evolution 
of isolat ed planetary sys tems. To this end, we use the MER- 
CURY (|Chamberd|l9 99) package to model the evolution of 
the system in between two encounters. MERCURY is a hy- 
brid integrator that allows symplectic integration of plane- 
tary systems over longer timescales, and can accurately deal 
with planet-planet scattering. The simulations typically con- 
serve energy at the \ AE/E\ = 1CT 6 - 10" 8 level. The switch 
from the CHAIN package to the MERCURY package and 
vice- versa is made when the encountering star is sufficiently 
far from the planetary system (r > 1000 AU), such that its 
effect on the outermost planet in the system is negligible. 

Star-planet and planet-planet collisions are modelled by 
merging two objects that enter each other's Roche lobe. In 
the case of a physical collision, they are replaced by a new 
particle with the combined mass of the two objects and with 
a position and velocity of their centre-of-mass. 

Throughout our simulations we continuously check for 
escaping planets, which are identified as follows. All planets 
at a distance r > -Riim from their host star are removed from 
the system, where i?n m is the average separation between the 
members of the star cluster. Planets at a distance 500 AU < 
r < -Riim are removed when (i) their binding energy E with 
respect to the central star is positive, (ii) a distance to the 
central star larger than 20 000 AU, and (iii) their velocity 
vector points away from the central star (f-v > 0). After 
escaper removal we calculate the velocity at infinity for each 
escaper correcting for the potential energy at the moment 
of removal. 

To improve statistics we carry out simulations for an 
ensemble of a hundred planetary systems. Each simulation 
starts out with an identical planetary system with randomly 
drawn initial orbital phases. The properties of the encoun- 
tering stars and the time scales between two subsequent en- 
counters drawn randomly from the distributions mentioned 
above. 



3 RESULTS 

3.1 Final planet configuration 

A planetary system experiences many encounters during the 
time which it spends in the star cluster. The cumulative 
effect of these encounters may result in the ejection of one or 
more of the planets, star-planet and planet-planet collisions, 
or a reconfiguration of the planetary system (see § 12.31 for 
details on the procedures involved). 

After simulating a hundred runs of target planetary sys- 
tems within their lifetime in the open cluster (~ 100 Myr), 
we obtain statistical results by combining the outcomes. Fig- 
ure[T]shows the fraction of planets that is still bound to their 
host star, as a function of their initial semi-major axis. As 
expected, planets with the smallest initial semi-major axis 
are more likely to survive. For model 1 a clear decline in the 
survival chances is seen with increasing initial semi-major 
axis. More than 50% of the planets with a < 10 AU have 
survived, while less than 30% of planets with a > 40 AU are 
bound to the system at the end of the simulations. The high- 
est retention rate is found for planets with a = 1 AU, with a 
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Figure 1. Fraction of survived planets as a function of their initial 
semi- major axis, for model 1 (blue) and model 2 (black). Inner 
planets have a higher probability of surviving after the stellar 
encounters, while outer planets can more easily escape the system. 



survival probability of more than 70%. For model 2 we see a 
similar, but steeper trend. A relatively large fraction of the 
Jupiters survive, while most of the outer planets escape the 
system. This steeper trend is a result of the mass spectrum 
of the planets in model 2: it is more difficult to eject the 
high-mass Jupiters, while it is easier to remove Uranus and 
Neptune from the system. 

The fraction of planets that is still bound after many 
encounters within the lifetime of its parent cluster is undis- 
putedly the most important parameter in assessing the effi- 
ciency of the combined effect of external perturbations and 
internal secular evolution in planetary systems. We there- 
fore examine the statistics of planet retention for all our 
data and plot the results in Figure [5] In the equal- mass sys- 
tem (model 1), 5% of the systems have lost all their plan- 
ets, while for 13% of the systems all planets remain bound 
despite the repeated stellar encounters. The peak in the dis- 
tribution gradually shifts to lower numbers over time. In 
model 2, on the other hand, 12% of the stars have lost all 
their planets, while only 1% of the systems retain all their 
planets. It is worth noting that almost half of the systems 
only has one planet (Jupiter) left at the end of the simu- 
lations. As explained above, this is caused by the unequal 
planetary masses in this configuration. 

Figure [3] shows the fraction of remaining bound plan- 
ets decreases smoothly with time. The equal-mass system 
(model 1) looses roughly one planet per 100 Myr, while the 
planet ejection rate is roughly two planets per 100 Myr for 
model 2. The latter set shows a decrease in the escape rate 
over time due to the reduced number of planets in the sys- 
tem. At the end of the simulations, each host stars has, on 
average, less than three planets remaining for model 1, and 
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Figure 2. The distribution of number of planets still bound to 
the host star at the end of the simulation, for model 1 (blue) and 
model 2 (black). At the end of the simulations, only 13% and 1% 
of the host stars have all their five planets in orbit, while 5% and 
12% of the host stars have lost all their planets, for model 1 and 
model 2, respectively. 
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Figure 3. The average number of planets bound to the central 
star as a function of time, for model 1 (blue) and model 2 (black). 



on average slightly more than one planet for model 2. The 
downward trend is still present at the end of the simula- 
tions, which simply indicates that additional planets will be 
removed from the system if it were to remain part of a star 
clusters beyond ~ 150 Myr. 

Distant stellar encounters perturb each of the planetary 
orbits to a certain degree. This generally result in a slight 
increase in the semi-major axis, the eccentricity, and the 
inclination. The cumulative effect of many encounters can 
substantially perturb the outer planetary orbits, which can 
in turn result in strong planet-planet interactions. A number 
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Figure 4. Eccentricity versus semi-major axis of the bound plan- 
ets at the end of the simulations, for model 1 (top) and model 2 
(bottom). The different colours in the top panel indicate the initial 
semi-major axis of 1 AU (red), 2.6 AU (yellow), 6.5 AU (green), 
16.6 AU (blue), and 42.3 AU (black). In the bottom panel the red, 
yellow, blue, and green colours represent Jupiter, Saturn, Uranus, 
and Neptune. At the end of the simulations, a substantial number 
of planets has obtained a higher eccentricity and a wider orbit. 



of planets obtain orbits that ultimately result in a physical 
collision with the central star. Also, some of the planets es- 
cape from the planetary systems and become a free-floating 
planet in the star cluster after perturbations. 

Figures f4] and [S] show the result of eccentricity and in- 
clination distribution as a function of semi-major axis after 
100 Myr. Note that the planets in model 2 (the four gas 
giants in the Solar system) lack low-eccentricity and low- 
inclination orbits, which is a consequence of their non-zero 
initial conditions. The chimney-like structures illustrate that 
it is much easier to change the eccentricity and inclination 
than the semi-major axis. The reason for this is that it is 
easier for neighbouring planets (and encountering stars) to 
exchange angular momentum, rather than orbital energy. It 
is worth mentioning that planets with an inclination larger 
than 90 degrees have retrograde orbits, and that all planets, 
including the innermost ones, have a chance of obtaining a 
retrograde orbit. 

Figure[6]shows the distribution of inclination angles ver- 
sus eccentricities of all remaining planets in the ensemble of 
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Figure 5. Inclination versus semi-major axis for the bound plan- 
ets at the end of the simulations, for model 1 (top) and model 2 
(bottom). Many planets have obtained both a higher inclina- 
tion due to stellar encounters, although the vast majority still 
have i < 40°. Note that some planets have obtained retrograde 
(i > 90°) orbits. 



models. We can see that, as expected, perturbed orbits both 
obtain higher eccentricities and higher inclinations. Note 
that there is no clear correlation with initial semi-major axis 
(indicated by the colour of the dots). Almost all highly in- 
clined and retrograde orbits have a large eccentricity, which 
may ultimately lead to Kozai cycles. 



3.2 Ejected planets 

The ejection velocity distribution of the escaping planets is 
shown in Figure [7] Note that there is no strong correlation 
between the ejection velocity and the initial semi-major axis. 
There is a weak correlation between the initial semi-major 
axis and the escape velocity. For the stars in model 1, the 
escaped innermost planets have a smaller escape velocity, 
which can be attributed to the deeper potential well which 
they have to escape from. The opposite result is seen for 
the Jupiters in model 2, which have a larger escape velocity 
than the other planets. This apparent contradiction can be 
explained by the mechanism at which these planets escape. 
In model 1, the innermost planets escape almost exclusively 
as a result of planet-planet scattering. In model 2, on the 
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Figure 6. Inclination versus eccentricity of the bound planets 
at the end of the simulations, for model 1 (top) and model 2 
(bottom). 



other hand, the Jupiters, which are substantially more mas- 
sive than the other planets, only escape as a result of direct 
interactions with encountering stars (see also § 13. 4[) . 

We can see that many of the survived planets have rel- 
atively large eccentricities and inclinations (Figure p4] and 
Figure [5]). These highly eccentric and highly inclined orbits 
in multi-planet system are unstable and might trigger esca- 
pers or collisions in the next several million years once their 
angular momentum deficit is high enough. 

From the range in time at which the planets escape, 
indicated with the horizontal lines in Figure we can 
see that the first escapers among the innermost planets 
in model 1 reach the boundaries of the system only after 
roughly 50 Myr. The escape velocity distribution among the 
other four planets is roughly constant in time and indepen- 
dent of their initial semi-major axes. This strongly suggests 
that escape among the innermost planets is primarily the re- 
sult of the secular evolution of the planetary system, rather 
than the direct result of the stellar encounters. For model 2 
we observe that each of the planets escape at random times, 
which indicates that at least a fraction of each of the planets 
is ejected due to direct interactions with encountering stars. 

The planets that are ejected from the system become 
individual star cluster members. Whether or not these free- 
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floating are immediately ejected from the star cluster as well, 
depends on their ejection velocity v e , the distance to the 
star cluster centre r. To a lesser extent, it also depends on 
the direction of the v e and the particle distribution within 
the cluster. Upon ejection from their planetary system, the 
planets obtain a velocity v = v e +v s , where v s is the velocity 
of the host star in the star cluster. The escape velocity w osc 
at a distance r from the centre of a star clusters is given by 



10.0000 



211 escapers in total 



2GM(r) 



(13) 



where G is the gravitational constant and M(r) is the en- 
closed mass within a radius r. 

When a free-floating planet remains bound to the 
star cluster, it may be re- captured by another stars 
l|Perets fe Kouwenhoven! I2012T ). or escape at a later time 
throu gh ejection or ev aporation. Using a microlensing sur- 
vey, ISumi et al.l (|201ll ) estimate that free-floating planets 
are almost twi ce as abundant as main se quence stars in the 
Galactic field. IVeras fc Raymond! (|2012h find that planet- 
planet scattering cannot explain the free-floating planet pop- 
ulation inferred from microlensing studies, and that other 
explanations, such as stripping of planets in star clusters, 
and stellar evolution, may well be necessary to explain the 
free-floating planet population. 



3.3 Star-planet and planet-planet collisions 

The perturbations induced by the encounters can result in 
instabilities in the planetary system which may result in 
star-planet and planet-planet collisions. We list the frac- 
tion of planets that experience collisions in Table [3] The 
error 8 on these fractions is 8 ~ N~ 1 d\/d~ 1 + AT -1 , where 
d is the number of ejections, collisions, or survivors, and 
N = 100 is the total n umber of objects in the sample 
(e.g.. lKouwenhoven|[2006l , p. 108). Here we have made the 
large-number assumption, and find that the error typically 
amounts to 8 — 5 — 10% for the data listed in the table. 

Table|3]shows that in the equal-mass case (model 1), the 
inner planets have a substantially larger probability of ex- 
periencing a physical collision than the outer planets, which 
are more frequently ejected from the system. For the inner- 
most planet, the number of escapers (15%) and the number 
of collisions with the central star (12%) are comparable. The 
vast majority of the collisions are between planets and the 
central star, while the collisions between two planets are 
rarely found in our simulations. The only case of a planet- 
planet collision is caused by a rare close encounter. Physical 
collisions are less frequent for model 2, where only 3% of 
the Saturns and 3% of the Neptunes collide with the central 
star. Collisions between Jupiter and the central star do not 
occur, as none of the other, lower-mass planets in the system 
are able to perturb Juptiter's orbit enough. 

Among additional tests with low-velocity encountering 
star in a near parabolic orbit with closest approach within 
100 AU, we find that planet-planet collision can occasionally 
occur. These results demonstrate that planet-planet inter- 
actions (rather than direct interactions with encountering 
stars) play an important role for the long-term evolution of 
perturbed planetary systems. 
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Figure 7. The velocities of the ejected planets as a function 
of time, for model 1 (top) and model 2 (bottom). Most of the 
inner planets (red data points) have rather low escape velocities, 
as compared to the high escape velocities of the outer planets 
(yellow, green, blue and black data points). The crosses indicate 
for each planet the time interval at which the planets escape and 
the standard deviation of their escape velocities. 



a initial 


Escapers 


Collisions 


Collisions 


Survivors 






w/ star 


w/ planets 




1.0 AU 


15 % 


12 % 


% 


73 % 


2.6 AU 


31 % 


7 % 


% 


62 % 


6.5 AU 


43 % 


5 % 


% 


52 % 


16.6 AU 


51 % 


2 % 


% 


47 % 


42.3 AU 


71 % 


1 % 


% 


28 % 


5.2 AU 


15 % 


% 


% 


85 % 


9.6 AU 


56 % 


3 % 


% 


41 % 


19.2 AU 


91 % 


% 


% 


9 % 


30.1 AU 


88 % 


3 % 


% 


9 % 



Table 3. The outcome of the multi-planet simulations for model 1 
(top) and model 2 (bottom). Note that that inner planets have 
substantially smaller survival chances than in the single-planet 
case (see Table U). 
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^initial 


Esc&pcrs 


Collisions 


Collisions 


Survivors 






w / stfir 


w / planets 




1.0 AU 


3 % 


% 


N/A 


97 % 


2.6 AU 


7 % 


% 


N/A 


93 % 


6.5 AU 


15 % 


% 


N/A 


85 % 


16.6 AU 


46 % 


% 


N/A 


54 % 


42.3 AU 


77 % 


1 % 


N/A 


22 % 


5.2 AU 


14 % 


% 


N/A 


86 % 


9.6 AU 


24 % 


% 


N/A 


76 % 


19.2 AU 


53 % 


% 


N/A 


47 % 


30.1 AU 


60 % 


% 


N/A 


40 % 



Table 4. The evolution of single-planet systems in an open cluster 
environment for model 1 (top) and model 2 (bottom). Note the 
collisions here exclusively represent mergers between a planet and 
the central star. 

3.4 Multi-planet systems versus single-planet 

systems 

We study the evolution of multi-planet systems in dense 
stellar environments. In order to demonstrate the impor- 
tance of multiplicity and planet-planet scattering in these 
systems, we carry out additional simulations in which we 
follow the evolution of single-planet systems in identical en- 
vironments. In this experiment we carry out simulations in 
which the planetary orbits are as described in § 12.21 but in 
this case only with for stars with a single planet. We use 
exactly the same properties of the encountering stars and 
encounter intervals, such that we are able to make a proper 
comparison. 

The results for the evolution of the single-planet sys- 
tems are listed in Table [4] The differences with those of the 
multi-planet systems (Table [3| are obvious: planet-planet 
scattering is of great importance to the evolution of multi- 
planet systems. Among the close-in planets in model 1 we 
observe many more escapers and collisions in the multi- 
planet case, with respect to the single-planet case. The 
differences are smaller for the outermost planets, as these 
are mostly affected by direct interactions with encounter- 
ing stars. Our data hints that outer planets have a slightly 
higher survival chance in multi-planet systems. One possible 
explanation for this is the interaction with the inner plane- 
tary system, which may exchange angular momentum with 
a strongly perturbed outer planet, which may mildly reduce 
its chances of escape. However, statistics should be improved 
in order to make a clear statement about this difference. In 
the case of the four Solar system planets (model 2), a similar 
trend is observed. One important difference is that, although 
the evolution of the outer three planets is substantially af- 
fected by the multiplicity, the Jupiters are unaffected. The 
reason for this is that the Jupiters are much more massive 
than the other planets, and therefore relatively unaffected 
by planet-planet scattering. For this reason, the surviving 
fraction of Jupiters is roughly 85%, irrespective of whether 
other, lower-mass planets are present in the system. 



4 CONCLUSIONS 

We have studied the evolution of planetary systems in star 
clusters similar to the ONC over an evolutionary time of 



the order of 10 years. We study two types of planetary sys- 
tems in particular. In the first set of simulations (model 1) 
we study the effect of encounters between planetary systems 
consisting of five Jupiter-mass planets in the semi-major axis 
range of 1 — 50 AU and single stars that approach these 
planetary systems within a distance less than 1000 AU. In 
the second set of simulations, we consider systems similar 
to our own Solar system, containing the Sun, Jupiter, Sat- 
urn, Uranus, and Neptune. The interactions are integrated 
using the classical Aarseth-Mikkola chain integrator in reg- 
ularised coordinates with high efficiency and accuracy. For 
the time interval between two subsequent interactions the 
MERCURY integration package is used to follow the inter- 
nal secular evolution of the planetary systems. Our main 
results are as follows: 

(i) In multi-planet systems in open star clusters, planets 
of all periods are either directly or indirectly affected by 
stellar encounters. The inner planetary orbits are affected 
as a result of interaction with perturbed outer planets. 

(ii) For the equal-mass system (model 1), the planets 
starting at a w 42 AU, approximately 25% survives, while 
~ 75% escapes, both in the single-planet and multi-planet 
cases. For those planets starting out at a = 1 AU in a 
multi-planet system, approximately 15% escapes, 12% col- 
lides with the central star, and 73% survives. In the single- 
planet case, however, 93% of the planets at a = 1 AU sur- 
vives, 3% escapes, and we observe no collisions. 

(iii) For the four giant planets in our Solar system 
(model 2), we observe a similar trend, although we see that 
Jupiter is not substantially affected by planet-planet en- 
counters, because of its large mass with respect to the other 
three planets. 

(iv) The planetary systems loose their planets at a rate of 
one per 100 Myr (model 1) and two per 100 Myr (model 2), 
respectively. A gradual decrease in the escape rate is seen 
over time, due to the smaller number of remaining planets 
in each system. 

(v) Although our results apply to the two specific plane- 
tary systems a given star cluster environment, these results 
do demonstrate that multiplicity cannot be ignored when 
studying the evolution of planetary systems in crowded en- 
vironments. 

(vi) Most of the planets' semi-major axes are increased 
slightly after the each encounters. The majority of the plan- 
ets obtain inclined and eccentric orbits. Generally, these ef- 
fects are strongest for the outer planets. Typically a few 
percent of the planets are found in highly inclined or retro- 
grade orbits at the end of the simulations. This applies to 
planets of any initial semi-major axis, and can be attributed 
to both the stellar encounters and planet-planet scattering. 
The latter process is dominant for the inner planets. 

(vii) The escaping planets are ejected with velocities up 
to several kms -1 , but, as Figure [7] shows, the majority of 
these newly formed free-floating planets are unable to obtain 
velocities beyond the escape velocity of their host star clus- 
ter. In addition, we find that the escape velocity increases 
with the initial semi-major axis of the planets. 

In summary, the presence of multiple planets in a plane- 
tary system results in a significantly different evolution as 
compared to single-planet systems, as a result of planet- 
planet interactions. Our results on the orbital evolution and 
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also our planet liberation rates per crossing time are in 
fair agreement with t he earlier direct cluster simulations of 
ISpur zcm et all l|2009h . who find in their Table 3 approxi- 
mately one planet liberation per crossing time for their wide 
planets (up to 50 AU). 



5 DISCUSSION 

Our model is a simplification of reality in a number of ways, 
and our assumptions should be improved and addressed in 
future work, ideally through directly modelling multi-planet 
systems in star clusters. 

We draw masses from an invariant Chabrier mass dis- 
tribution (Eq.rjJ), we draw the relative velocities from a local 
Maxwellian distribution function, based on the average ve- 
locity dispersion of the Plummer model (Eq.(3|, we draw im- 
pact parameters from a linear probability distribution func- 
tion (Eq.[5]), and we assume that subsequent encounters are 
uncorrelated, well separated in time (Eg. Hip. In a real evolv- 
ing star cluster all these assumptions are to some degree ap- 
proximations, but the full star cluster simulation is a much 
larger numerical effort than our work here, and to some ex- 
tent it is more complex to disentangle the effects of individ- 
ual encounters in a full simulation, as can be seen from the 
rather complicated procedures used in ISpurzem et~all l|200gft 
to identify encounters in TV-body simulations. It should be 
noted, however, that mass segregation, core collapse and the 
finite size of the core of a star cluster will lead to a time vari- 
ation of the encounter parameters, and to a deviation from 
the simple distribution functions assumed here. In particu- 
lar, it is known that in real star clusters the ve locity dis- 
tribution function can be far from a Maxwellian (|Plummerl 
ll91ll : ICohrJll98Ch . 

The encounter rate in a star cluster decreases with time 
as a result of expansion. In addition, a star may escape from 
a star cluster within a short time, after which the effect 
of stellar encounters on their planetary systems are negli- 
gible. We obtai n approximate evolutionar y tracks using the 
EMACSS code (I Alexander fc Gielesl[2012l ). and estimate the 
stellar encounter rate as a function of time under the ap- 
proximation that the interactions are domin ated by gravi- 
tational focusing (see iMalmberg et al.ll20l"lh . As expected, 
the encounter rate gradually decreases over time, by a factor 
of two after 30 Myr down to a factor of six after 100 Myr. 
These estimates demonstrate that we may have overesti- 
mated the number of encounters typically experienced by a 
planetary system. On the other hand, star clusters may form 
sub-virial, and star forming regions observed to form with 
a significant amount of substructure in both the position 
and the velocity distribution, which can boost the number 
of close e ncounters during th e early stages of star cluster 
evolution j Allison et al.l 12009). Finally, if a star cluster fol- 
lows the evolutionary track listed above, many of its stars 
will suffer additional encounters until the cluster is com- 
pletely dissolved. This shows again the need for full TV-body 
simulations, which naturally include all processes described 
above. Despite these approximations, we have demonstrated 
the importance of multiplicity of planetary systems, and we 
expect that improved future simulations will give qualita- 
tively the same result. 

We have also neglected the presence of binary 



systems, despite their p resen ce in star clusters (e.g., 
iKouwenhoven et al.l 120051. l2007h and in the Galactic field 
l|Duquennov fe Mayor! Il99ll ). The effect of binary systems 
is complex, and depends strongly on the semi-major axis 
a of the binary and the distance of closest approach p of 
the centre-of-mass of the binary. Encounters with sufficiently 
tight binaries (a <C p) have a very similar effect as encoun- 
tering single stars, since they can, to first order, be consid- 
ered as single bodies with a larger mass (up to a factor two 
for equal-mass binary systems). Very wide binaries (a ^> p) 
also have a similar effect as single stars, as a planetary sys- 
tem only feels the nearby component of the binary system. 
The main difference with the encounters of single stars is 
that, while the centre-of-mass of a wide binary approaches 
impact parameter b, the periastron distance of the nearest 
binary component may come closer to the planetary system. 
Complex situations arise in the case of intermediate-period 
binaries with a as p. In this case, both components of the 
binary system can affect the planetary system, but more 
importantly, highly destructive three-body interactions can 
occur. Depending on the properties of the binary popula- 
tion in a star clusters, these types of interactions may or 
may not have an important effect on the survival of plane- 
tary systems. Further work is definitely necessary to study 
the effect of binary systems. However, from the arguments 
above one can see that the inclusion of binaries, whether 
they are tight, intermediate, or wide, will generally cause 
further destruction of planetary systems. 

Finally, it would be interesting to see how our results 
change when considering the immense diversity among the 
observed exoplanet systems. Systems with a different num- 
ber of planets, planetary masses, and relative separations 
will result in different outcomes. However, the general out- 
come of our simulations (planets at any orbital separation 
are either directly or indirectly affected by distant stellar 
encounters) is unlikely to change. 

Future work will combine our ansatz with full star clus- 
ter simulations. The secular evolution of planetary systems 
between encounters in an TV-body simulation of a large star 
cluster can be treated by either an adapted version of MER- 
CURY code or by a generalisation of Lagrange's planetary 
equations, while encounters will take place as in a real star 
cluster initialised by the time-dependent environment of a 
full TV-body simulation. 
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